function [obj, Delta_xi, S, v_c, v_m, s1_predicted, s1_predicted_y, pr_stay,...
    s1_predicted_y_c, s1_predicted_y_m,...
    moment_mean_all, naic_xi,...
    var_cov, N_ist_unmerge] ...
    = demand_objective_inner(x, K, D, N,...
    pr_fc, x_fc, lapse, beta_c, crra, rho, sigma_c1, sigma_c2, eulerc,...
    N_y, pr_y, y_medi, T1, T2, B1_c, B2_c, naic_unique,...
    no_insurer_fe, ind_select, delta, use_log_n, myopic,...
    resource_y1, resource_y2)
%% Preliminaries
% demand parameters to be estimated
alpha      = x(1);     
lapse_u    = x(2);     
variety_u  = 1;     
xi         = x(4:end); 

% inlude insurer's naic_code
naic_xi = [naic_unique xi];

% get data on prices
p1 = D.p1; % N.ist x 1
r2 = D.r2; % N.ist x K

% get data on 1st period market shares
s1_data = D.s1; 
ind_f = find(D.major==0);
s1_data(ind_f) = D.s1(ind_f)./D.n_insurers(ind_f);

% get data on # of plans/insurers in each product
n_variety = D.n_insurers;

% tolerance level for contraction
tol=0.001;

% initial guess of Delta_xi
if isfile('Delta_xi_interim.mat')==1
    temp = load('Delta_xi_interim.mat');
    Delta_xi = temp.Delta_xi;
else
    Delta_xi = zeros(N.ist,1);
end

%% Consumer's value from not purchasing any insurance
% # of states x (1+K). 1st col has the numeric state value
state_pr_s = unique([D.state D.pr_s], 'rows');

% loop thru states
for ss = 1:N.s
    
    % pr_s: K x 1
    pr_s = state_pr_s(ss,2:K+1)';
    
    [v0, U0] = consumer_v0_ftn(N_y, rho, crra, beta_c, K, ...
        pr_s, x_fc, pr_fc, y_medi, B1_c, B2_c, T1, alpha,...
        lapse_u, resource_y1, resource_y2);
    
    S(ss).v0 = v0;
    S(ss).U0 = U0;
end

%% Solve the contraction mapping
eps = 100;
iter = 1;

while (eps>tol)
    
    %% Consumer's choice-specific value from purchasing a contract    
    correct_beliefs=1;
    [v_c, pr_stay] = consumer_v_ftn(p1, r2, N, naic_xi, lapse, N_y, ...
        rho, crra, beta_c, K, state_pr_s, S, sigma_c2, eulerc, B1_c, B2_c, T1, ...
        Delta_xi, alpha, D,...
        no_insurer_fe, ...
        n_variety, variety_u, delta, use_log_n, correct_beliefs,...
        resource_y1, resource_y2);
    
    correct_beliefs=0;
    [v_m, ~] = consumer_v_ftn(p1, r2, N, naic_xi, lapse, N_y, ...
        rho, crra, beta_c, K, state_pr_s, S, sigma_c2, eulerc, B1_c, B2_c, T1, ...
        Delta_xi, alpha, D,...
        no_insurer_fe, ...
        n_variety, variety_u, delta, use_log_n, correct_beliefs,...
        resource_y1, resource_y2);

    %% Derive predicted market shares
    % store each insurer's predicted share in a given market
    s1_predicted = zeros(N.ist,1);
    
    % store each insurer's predicted share in a given market from each
    % consumer type -> needed for supply side estimation only
    s1_predicted_y = zeros(N_y, N.ist);
    s1_predicted_y_c = zeros(N_y, N.ist); % when correct beliefs
    s1_predicted_y_m = zeros(N_y, N.ist); % when misbeliefs 
    
    for mm = 1:N.st
        % unmerge v_c and v_m
        ind_m = find(D.mkt==mm & D.major==1)'; % 1,...,N.ist
        ind_f = find(D.mkt==mm & D.major==0)'; % 1,...,N.ist

        n_fringe = D.n_insurers(ind_f);
        n_major = length(ind_m);
       
        v_c_m = v_c(:, ind_m);
        v_c_f = v_c(:, ind_f);
        v_c_unmerged = [v_c_m repmat(v_c_f, 1, n_fringe)];

        v_m_m = v_m(:, ind_m);
        v_m_f = v_m(:, ind_f);
        v_m_unmerged = [v_m_m repmat(v_m_f, 1, n_fringe)];

        % market mm's state
        state_mm = unique(D.state(ind_m));
        
        % consumers' value from staying uninsured
        v0_mm = S(state_mm).v0;  % N_y x 1
        
        NN_unmerged = size(v_c_unmerged,2);

        % choice probabilities conditional on using actual prices 
        ind_temp = ones(1,NN_unmerged);
        s1_ind_c = share_ind_ftn(v_c_unmerged, ind_temp, v0_mm, sigma_c1); %N_y x (N_e+1)
        
        % choice probabilities conditional on using initial rates for future rates
        s1_ind_m = share_ind_ftn(v_m_unmerged, ind_temp, v0_mm, sigma_c1); %N_y x (N_e+1)
        
        % N_y x (N_e+1)
        s1_ind = (1-myopic)*s1_ind_c + myopic*s1_ind_m;

        % probability of purchasing each of N_e products
        s1_ind1 = s1_ind(:,1:end-1); % N_y x N_e

        s1_predicted_y(:, ind_m) = s1_ind1(:,1:n_major); % N_y x n_major
        s1_predicted_y(:, ind_f) = s1_ind1(:,n_major+1); % N_y x 1

        % 1st period mkt share,
        s1_mkt = sum(s1_ind1.*repmat(pr_y, 1, NN_unmerged))'; % N_e x 1

        s1_predicted(ind_m') = s1_mkt(1:n_major);
        s1_predicted(ind_f') = s1_mkt(n_major+1);
        
        % store individual purchase prob conditional on correct & misbeliefs
        s1_ind1_c = s1_ind_c(:,1:end-1); % N_y x N_e
        s1_predicted_y_c(:, ind_m) = s1_ind1_c(:,1:n_major);
        s1_predicted_y_c(:, ind_f) = s1_ind1_c(:,n_major+1);

        s1_ind1_m = s1_ind_m(:,1:end-1); % N_y x N_e
        s1_predicted_y_m(:, ind_m) = s1_ind1_m(:,1:n_major);
        s1_predicted_y_m(:, ind_f) = s1_ind1_m(:,n_major+1);
        
    end
    
    %% Contraction mapping
    % compute error term
    eps = max(abs((s1_data./s1_predicted) - 1));
    
    % contraction mapping
    if no_insurer_fe==0
        log_share_diff = (log(s1_data) - log(s1_predicted))/B1_c;
        Delta_xi_new   = Delta_xi + log_share_diff;

        % update
        Delta_xi = Delta_xi_new;
        
    elseif no_insurer_fe==1
        log_share_diff = (log(s1_data) - log(s1_predicted))/(B1_c + (beta_c^T1)*B2_c);
        Delta_xi_new   = Delta_xi + log_share_diff;

        % update
        Delta_xi = Delta_xi_new;
    end
    
    iter = iter+1;
    
end

%% Save Delta_xi to be called in the next function evaluation
save('Delta_xi_interim.mat', 'Delta_xi')

%% Delta_xi: deduct insurer FE and time FE
% compute insurer FE
xi_bar = zeros(N.ist,1);

for ii=1:length(naic_unique) 
    naic_ii = naic_unique(ii);
    ind = find(D.naic==naic_ii);
    if ii==1
        xii1 = mean(Delta_xi(ind));
    else
        xi_bar(ind) = mean(Delta_xi(ind)) - xii1;
    end
    
end

% compute time FE
xt_bar = zeros(N.ist,1);
year_unique = unique(D.year);

for tt=1:length(year_unique)
    year_tt = year_unique(tt);
    ind = find(D.year==year_tt);
    if tt==1 
        xt1 = mean(Delta_xi(ind));
    else 
        xt_bar(ind) = mean(Delta_xi(ind)) - xt1;
    end
    
end

% substract insurer FE and time FE
Delta_xi_demean = Delta_xi - xi_bar - xt_bar;

%% Moment conditions
% unmerge fringes: rows change from N.ist to NN:= (N.ist-N_st) + sum(D_nF)
N_st = N.st;
D_major = D.major;
if sum(D_major==0) ~= N_st
    disp('ERROR: sum(D_major==0) ~= N_st')
end
D_nF = D.n_insurers(D_major==0); % N_st x 1, contains # of fringes in a given mkt

arg_share=0;
D_unmerge_Delta_xi_demean = unmerge_fringe(Delta_xi_demean, D_major, arg_share, D_nF, N_st); % NN x 1
D_unmerge_iv_prem_nstate = unmerge_fringe(D.iv_prem_nstate, D_major, arg_share, D_nF, N_st); % NN x 1
D_unmerge_delta_approved_cum = unmerge_fringe(D.delta_approved_cum, D_major, arg_share, D_nF, N_st); % NN x 1

% sample size
N_ist_unmerge = length(find(isnan(D_unmerge_iv_prem_nstate)==0));

% moments
ind = find(isnan(D_unmerge_iv_prem_nstate)==0);
mmt_prem_nstate = D_unmerge_Delta_xi_demean(ind).*(D_unmerge_iv_prem_nstate(ind)-mean(D_unmerge_iv_prem_nstate(ind)));

ind = find(isnan(D_unmerge_delta_approved_cum)==0);
mmt_delta_increase_reg = D_unmerge_Delta_xi_demean(ind).*(D_unmerge_delta_approved_cum(ind)-mean(D_unmerge_delta_approved_cum(ind)));

% choose moments to use
MMT(1).vec = mmt_prem_nstate;
MMT(2).vec = mmt_delta_increase_reg;
N_moments = 2;

% means
moment_mean_all = zeros(N_moments,1);
for nn=1:N_moments
    moment_mean_all(nn) = sum(MMT(nn).vec)/N_ist_unmerge;
end

% variance-covariance matrix: diagonal entries
var_cov = zeros(N_moments, N_moments);
for nn=1:N_moments
    hh = MMT(nn).vec';
    var_cov(nn,nn) = (1/N_ist_unmerge)*(hh*hh');
end

% variance-covariance matrix: off diagonal entries
ind = find(isnan(D_unmerge_iv_prem_nstate)==0 & isnan(D_unmerge_delta_approved_cum)==0);
hh = (D_unmerge_Delta_xi_demean(ind).*(D_unmerge_iv_prem_nstate(ind)-mean(D_unmerge_iv_prem_nstate(ind))))...
    .*(D_unmerge_Delta_xi_demean(ind).*(D_unmerge_delta_approved_cum(ind)-mean(D_unmerge_delta_approved_cum(ind))));
var_cov(1,2) = (1/N_ist_unmerge)*sum(hh);
var_cov(2,1) = var_cov(1,2);

% objective function
obj = moment_mean_all'*(var_cov\moment_mean_all); % optimal weight matrix

%% For fringes' market share, store agg values
ind_f = find(D.major==0);
s1_predicted(ind_f) = s1_predicted(ind_f).*D.n_insurers(ind_f);
s1_predicted_y(:,ind_f) = s1_predicted_y(:,ind_f).*repmat(D.n_insurers(ind_f)', N_y, 1);
s1_predicted_y_c(:,ind_f) = s1_predicted_y_c(:,ind_f).*repmat(D.n_insurers(ind_f)', N_y, 1);
s1_predicted_y_m(:,ind_f) = s1_predicted_y_m(:,ind_f).*repmat(D.n_insurers(ind_f)', N_y, 1);

